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We consider a Hubbard model on a square lattice with an additional interaction, W, which depends 
upon the square of a near-neighbor hopping. At half-filling and a constant value of the Hubbard 
repulsion, increasing the strength of the interaction W drives the system from an antiferromagnetic 
Mott insulator to a d x 2_ y 2 superconductor. This conclusion is reached on the basis of zero temper- 
ature quantum Monte Carlo simulations on lattice sizes up to 16 x 16. 
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At half-filling, the two-dimensional (2D) Hubbard 
model on a square lattice with a near-neighbor one- 
electron hopping t and an onsite Coulomb interaction 
U has an insulating ground state with long range anti- 
ferromagnetic (AF) order. When carriers are doped into 
this insulating state, the possibility of a transition to a 
d x 2_ y 2 superconducting state has been extensively pur- 
sued. However, at energy scales and lattice sizes accessi- 
ble to numerical simulations of the 2D Hubbard model, 
no definite sign of superconductivity has been found in 
the vicinity of the Mott transition yj. An alternate ap- 
proach for studying the transition from the Mott insu- 
lator to a d x 2_ y 2 superconducting state is to keep the 
system at half-filling, but add an interaction which can 
drive the system into a d x 2_ y 2 superconducting ground 
state. Here, we introduce a simple form for such an in- 
teraction which depends on the square of a near-neighbor 
hopping. We note that it can formally be obtained from 
a Su-Schrieffer-Heeger Q interaction in the antiadiabatic 
limit. It may also be viewed as derived from purely elec- 
tronic origins near the AF Mott insulator. However, our 
purpose is not aimed at describing the pairing mechanism 
but rather at studying the transition between an antifer- 
romagnetic Mott insulating state and a d x 2 _ y 2 supercon- 
ducting state. In this sense, we view the interaction as 
a formal way of obtaining a model which exhibits this 
transition and is suitable for Monte-Carlo simulations in 
that there are no fermion sign problem. 

The basic half-filled Hubbard model that we will study 
has the Hamiltonian 
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Here, cl { c % a ) creates (annihilates) an electron with 



component of spin a on site i, n-? 
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energy will be measured in units of t. 
that we will add has the form: 
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The interaction 
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with positive W. We note that the interaction Hw can 
be generated from a Su-Schrieffer-Heeger Q term with 
Einstein oscillators: 
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Integrating out the phonons and taking the antiadiabatic 
limit (M -> 0), generates H w with W = A t L>" 1 A/2. 
Pairing mechanism along those lines were considered in 
HQ. This derivation should be compared to the case of 
the attractive Hubbard model which may be obtained in 
a similar way provided that electron-phonon interaction 
is chosen to be of the Holstein form ||. In fact, one of 
the reasons for choosing this form is that it has a simple 
Hubbard-Stratonovitch representation which is useful in 
constructing the Monte-Carlo simulation. 
The Hamiltonian 



H = H r 
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has the possibility of exhibiting a quantum transition be- 
tween an AF Mott insulating state and a superconduct- 
ing phase. When W — 0, the half filled Hubbard model 
with a finite U is known to be a Mott insulator with long- 
range AF order. The interaction H\y can be decomposed 
into single-particle next-nearest neighbor hopping terms, 
singlet pair-hopping terms involving on-site singlet pairs, 
triplet pair-hopping processes, and terms of the form: 
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where Ai 



(ct cl ~ 



ctcl p ) /y/2. In the pres- 
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ence of the Hubbard Hamiltonian Hjj, the latter term 
(||) dominates the low energy physics, and as we will see, 
leads to a d x 2_ y 2 superconducting state. 

We have carried out our simulations with a zero tem- 
perature quantum Monte-Carlo (QMC) algorithm ^,0. 
At half- filling the Hamiltonian, Eq.(||), is particle- hole 
symmetric at any U and W so that there is no sign prob- 
lem and ground state simulations on large lattices (up to 
16 x 16) were carried out without any complications. This 
statement is valid even for twisted boundary conditions 
for which 
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with $o = hc/e the flux quanta and L the linear length 
of the square lattice. The boundary conditions given by 
Eq. (Q) account for a magnetic flux threading a torus on 
which the lattice is wrapped. To take advantage of the 
efficiency of a single spin-flip algorithm, we have carried 
out a discrete Hubbard Stratonovitch transformation of 
the Hp term. The transformation produces an overall 
systematic error of order (AtW) 3 where At corresponds 
to the imaginary time step entering the path integral 
formulation. This systematic error is negligible compared 
to the one produced by the Trotter decomposition which 
is of order (At) 2 . Unless mentioned otherwise, we have 
carried out our simulations at Art — 0.0625. The details 
of the algorithm will be presented elsewhere Q . 

We first concentrate on the charge degrees freedom 
with U/t = 4, and study the ground state as a func- 
tion of W/t. In order to do so, we have computed the 
ground state energy as a function of the twist in the 
boundary condition: Eq((&). For an insulator, the wave 
function is localized and hence, an exponential decay of 
AE (&) = E (&) - E (& /2) as a function of lattice size 
is expected Q. In the spin density wave (SDW) approx- 
imation for the half-filled Hubbard model, one obtains 
AEq($) = a(&)L exp (— L/£) where £ is the localization 
lenght of the wavefunction. On the other hand, for a 
superconductor, AEq(<&) shows anomalous flux quantiza- 
tion: AEq ($) is a periodic function of $ with period $o/2 
and a non vanishing energy barrier is to be found between 
the flux minima [jTo|— p~2|| so that Ai? ( < f ) o/4) remains fi- 
nite as L — * oo. Fig. la shows AEo(Qq/4) versus 1/L for 
various values of W/t. One observes a change in the size- 
scaling of A£'o( < i ) o/4) as W/t decreases from W/t = 0.5 
to W/t = 0.22. From these measurments, we estimate 
that the change occurs in the vicinity of W/t = 0.3. For 
values of W/t < 0.3 AE ($o/4) is consistent with the 
SDW form whereas for W/t > 0.33 AE ($ /4) may be 
fitted to a 1/L form and scales to a finite value. The 
extrapolated value of AEq(&o/4) versus W/t is plotted 
in Fig. lb and the quantum transition between a Mott 
insulator and superconductor occurs at W/t ~ 0.3. 

In order to determine the symmetry of the order pa- 



rameter in the superconducting state, we have calculated 
pair- field correlations in the s and d x 2_ y 2 channels: 



with 
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Here, f s (S) = 1 and f d (S) = 1(-1) for S = ±a x (±a y ). 
Fig. 2 shows plots of P d . s (L/2, L/2) for W/t = 0.6 where 
the system is superconducting and W/t = 0.1 where it is 
not. At W/t — 0.6, the dominant signal at long distances 
(L = 16) is obtained in the d x 2_ y 2 channel. 

At the mean-field level, the symmetry of the order pa- 
rameter will determine the functional form of the single 
particle occupation number, n(k). For a d x 2_ y 2 super- 
conductor the BCS result yields: 



n{k) 
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where = — 2t (cos(k x ) + cos(k y )) and Ar — 

A (cos(k x ) — cos(k y )) . As apparent from Eq. (|10| ) in 
the k = A;(l,l) direction the d x 2_ y 2 gap vanishes and 
n(k) shows a jump at the Fermi energy, whereas in the 
k = k(l, 0) direction n{k) is a smooth function of k. Pre- 
cisely this behavior in n(k) may be detected in the QMC 
data at W/t = 0.6 as shown in Fig. 3a. For comparison, 
we have plotted n(k) at W — where it is expected to 
scale to a smooth function in the thermodynamic limit 
(see Fig.3b). 

We now consider the spin degrees of freedom. We have 
computed the real space spin-spin correlations: 



S(f) = -(S(r)S(6)) 



(11) 



where S (r) is the spin operator at site r. For values 
of W/t < 0.3 and lattice sizes ranging from L = 4 
to L = 12, S(L/2,L/2), may be fitted to a 1/L form 
and scales to a finite value, as shown in Fig. 4a. We 
therefore conclude that long range AF order is present 
for W/t < 0.3. The associated staggered moment, 
m = liniL^oo y/3S(L/2,L/2), is plotted in Fig 4b. The 
data is consistent with a continuous decay of m as W/t 
increases towards 0.3. At W/t = 0.3, we were unable 
to distinguish m from zero within our statistical uncer- 
tainty. Hence, we conclude that long-range AF order 
vanishes at W/t ~ 0.3. Therefore, within our numeri- 
cal resolution, the antiferromagnetic transition point is 
not separated from the superconductor-insulator tran- 
sition point. Well within the d x 2_ y 2 superconducting 
phase the spin-spin correlations remain sizable. In fact, 
at W/t = 0.6, lattice sizes ranging from L — 4 to 
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L = 16, S(L/2, L/2) scales as L~~ a with a = 1.16 ± 0.01 
as shown in the inset of Fig. 4a. This slow decay of 
the spin-spin correlations in the superconducting state 
arise because of the nodes in the d x 2_ y 2 gap. Follow- 
ing Ref. 113], one can approximate the spin susceptibil- 
ity, x{Qi iu m ) , in the superconducting state by inserting 
the irreducible BCS spin susceptibility, xo(Qy i^m), in the 
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corresponds to Matsubara frequencies. Within this ap- 
proximation and at half-band filling, the spin-spin cor- 
relations for a d x 2_ y 2 superconducting order parameter 
show a powerlaw decay. In contrast, for an s— wave order 
parameter, an exponential decay is obtained. 

We have found a simple model, Eq. (||) , which exhibits 
a ground state transition from an AF Mott insulator to 
a d x 2_ y 2 superconductor. Specifically, for U/t — 4 at 
half- filling ((n) = 1) we have found that this transition 
occurs at W c ~ 0.3i. As a natural consequence of these 
results, we expect for values of W < W c the occurrence 
of a quantum transition between an AF Mott insulator 
and a d x 2_ y 2 superconductor as a function of chemical 
potential, \x. It has recently been concluded that for the 
2D Hubbard model (W = 0), the quantum transition be- 
tween the AF Mott insulator and metallic state driven 
by the chemical potential, belongs to a universality class 
characterized by an unusually large dynamical exponent, 
z = 4 [p^-^6[ . The presence of the Hw term may by-pass 
this suppression of coherence observed at W = and in- 
troduce another independent energy scale to drive the 
Mott transition. At constant U and in the W, /i plane, 
the Hamiltonian (j^) is expected to lead to a rich struc- 
ture of quantum transitions involving the Mott insulat- 
ing state, the d-wave superconducting state and metallic 
states as well as crossovers at finite temperatures. Al- 
though the justification of the model from a microscopic 
point of view is at present uncertain, it provides us with 
a model to study some of the salient ground state and 
finite temperature features of a d x 2_ y 2 superconductor 
at an energy scale accessible to numerical simulations. It 
also allows us to study a quantum transition between an 
AF Mott insulator and a d x 2_ y 2 superconductor. 
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Figure captions 

Fig. 1 (a) A£ ($o/4) = £o($o/4) - £ (*o/2) 
versus 1/L for several values of W/t. For 
W/t < 0.3 the solid lines correspond to a 
least square fit of the data to the SDW form: 
Lexp(-L/£). For W/t > 0.3 the QMC data 
is compatible with a 1/L scaling. The solid 
lines are least square fit to this form, (b) Ex- 
trapolated value of AE (3> /4) versus W/t. 
The solid line is a guide to the eye. 

Fig. 2 d x 2_ y 2 (triangles) and s-wave (circles) pair- 
field correlations versus 1/L. 

Fig. 3 (a) n{%) at W/t = 0.6, U/t = 4 and (n) = 1. 
Lattices form L = 8 to L — 16 were consid- 
ered, (b) same as (a) but for W/t = 0. The 
calculations in this figure were carried out at 
$ = (see Eq. g). 
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4 (a) S(L/2, L/2) versus 1/L for several values 
of W/t. The solid lines correspond to least 
square fits of the QMC data to the form 1/L. 
Inset: S(L/2,L/2) versus 1/L at W/t = 0.6. 
The solid is a least square fit to the form L~ a . 
(b) Staggered moment as obtained from (a) 
versus W/t. The data point at W/t — is 
taken from reference flTj]. At W/t = 0.3, we 
were unable to distinguish m from zero within 
our statistical uncertainty. The solid line is 
a guide to the eye. The calculations in this 
figure were carried out at $ = (see Eq. ^). 
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